import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

colors = [(0, 0, 0), (1, 0, 0)] # black to red
cm = LinearSegmentedColormap.from_list(
        "Custom", colors, N=256)

Gaussian Beams

The complex amplitude of Gaussian beam has the form,

$$ \begin{equation} U(r) = A_{0} \dfrac{W_{0}}{W(z)} exp\left[-\dfrac{\rho^2}{W^2(z)}\right] exp\left[-jkz - jk\dfrac{\rho^2}{2R(z)} + j\zeta(z)\right] \end{equation} $$

$ \textbf{Beam width} \; \begin{equation} W(z) = \sqrt{\dfrac{\lambda z_{0}}{\pi}} \sqrt{1 + \left(\dfrac{z}{z_{0}}\right)^2} = W_{0} \sqrt{1 + \left(\dfrac{z}{z_{0}}\right)^2} \\ \end{equation} $ $ \begin{equation} \textbf{Wavefront radius of curvature} \; R(z) = z \left[1 + \left(\dfrac{z_{0}}{z}\right)^2\right] \\ \end{equation} $ $ \begin{equation} \textbf{Gouy phase} \; \zeta(z) = tan^{-1}\left(\dfrac{z}{z_{0}}\right) \\ \end{equation} $ $ \begin{equation} \textbf{Rayleigh range} \; z_{0} = \dfrac{\pi W_{0}^2}{\lambda} \\ \end{equation} $ $ \begin{equation} \textbf{Beam waist} \; W_{0} = \sqrt{\dfrac{\lambda z_{0}}{\pi}} \end{equation} $

wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)


print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')


x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)

xx, yy = np.meshgrid(x, y)

rho = np.sqrt(xx**2 + yy**2)

Beam waist: 0.014192 mm
Divergence angle: 0.028385 rad
Depth of focus: 2.000000 mm

Beam width

Beam width at y-z plane

def W(z, wavelength, z0):
    W0 = np.sqrt(wavelength * z0 / np.pi)
    
    return W0 * np.sqrt(1 + (z / z0) ** 2)

_W = W(z, wavelength, z0)
_theta = np.sqrt(wavelength / (z0 * np.pi)) * z

fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(111)
ax.plot(z, _W, color='blue')
ax.plot(z, -_W, color='blue')
ax.plot(z, _theta, linestyle='--', color='red')
ax.plot(z, -_theta, linestyle='--', color='red')

ax.grid(True)
# ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$W(z)$", fontsize=15)
plt.show()

Radius of curvature

def R(z, z0):
    return z * (1 + (z0 / z) ** 2)

fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(111)
ax.plot(z, R(z, z0), color='blue')
ax.plot(z, z, linestyle='--', color='red')

ax.grid(True)
# ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$R(z)$", fontsize=15)
ax.set_ylim(-6*z0, 6*z0)
plt.show()

Intensity

The intensity of the Gaussian beam $I(\rho, z) = \left| U(r) \right|^2$ is,

$$ \begin{align} I(\rho, z) &= \left| A_{0} \right|^2 \left[\dfrac{W_{0}}{W(z)}\right]^2 exp\left[-\dfrac{2\rho^2}{W^2(z)}\right] \\ &= I_{0}\left[\dfrac{W_{0}}{W(z)}\right]^2 exp\left[-\dfrac{2\rho^2}{W^2(z)}\right] \end{align} $$

W0 = np.sqrt(wavelength * z0 / np.pi)
W0

0.014192480261642175

def intensity(rho, z, wavelength, z0, I0=1):
    W0 = np.sqrt(wavelength * z0 / np.pi)
    _W = W(z, wavelength, z0)
    
    return I0 * (W0 / _W) ** 2 * np.exp(-2 * (rho / _W) ** 2)
#     return I0 * (W0 / _W) ** 2 * np.exp(-2 * rho ** 2 / _W ** 2)

Intensity along beam axis $(\rho=0)$

I1 = intensity(0, z, wavelength, z0, I0=1)

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.plot(z, I1)
ax.grid(True)
ax.set_title(r"Intensity along beam axis ($\rho$=0)", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$I/I_{0}$", fontsize=15)
plt.show()

Intensity at the transverse plane

Assume z at 0, $z_{0}$, $2z_{0}$

z_list = [0, z0, 2*z0]

fig, ax = plt.subplots(1, 3, figsize=(24,5))

for i, _z in enumerate(z_list):
    I = intensity(rho, _z, wavelength, z0, I0=1)
    ax[i].imshow(I/I.max(), 
                 cmap=cm, 
                 extent=[x.min(), x.max(), y.min(), y.max()], aspect='equal')
    ax[i].set_title(f"z={_z:.3f} mm", fontsize=20)

Intensity at y-z plane

I = intensity(yy, z, wavelength, z0, I0=1)


fig = plt.figure(figsize=(8,2))
ax = fig.add_subplot(111)
ax.imshow(I, cmap='jet', origin='lower', extent=[z.min(), z.max(), y.min(), y.max()], aspect='auto')
# ax.imshow(I, cmap='jet', origin='lower')
ax.grid(True)
ax.set_title("Intensity at y-z plane", fontsize=20)
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel("y", fontsize=15)
plt.show()

def phase(z, rho, wavelength, z0):
    k = 2 * np.pi / wavelength
    longi = k * z - np.arctan(z / z0)
    trans = k * rho ** 2 / (2 * R(z, z0))
    
    return longi + trans


def amp(z, rho, wavelength, z0, A0=1):
    W0 = np.sqrt(wavelength * z0 / np.pi)
    _W = W(z, wavelength, z0)
    return A0 * (W0 / _W) * np.exp(-(rho / _W) ** 2)

z = np.linspace(-2*z0, 2*z0, 2000)
yy, zz = np.meshgrid(y, z)

U = amp(zz, yy, wavelength, z0, A0=1) * np.exp(-1j * phase(zz, yy, wavelength, z0))

fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
ax.imshow(np.abs(U.real).T**0.6, 
          cmap='jet', 
          origin='lower', 
          aspect='auto', 
          extent=[z.min(), z.max(), y.min(), y.max()])
ax.set_xlabel("z", fontsize=15)
ax.set_ylabel(r"$\rho$", fontsize=15)
plt.show()

Hermite-Gaussian beam

The complex amplitude of Hermite-Gaussian beam has the form,

$$ \begin{equation} U_{l,m}(x,y,z) = A_{l,m} \left[ \dfrac{W_{0}}{W(z)}\right] \mathbb{G}_{l}\left[\dfrac{\sqrt{2}x}{W(z)}\right] \mathbb{G}_{m}\left[\dfrac{\sqrt{2}y}{W(z)}\right] exp\left[-jkz - jk\dfrac{x^2+y^2}{2R(z)} + j(l+m+1)\zeta(z)\right] \end{equation} $$

$ \textbf{Hermite-Gaussian function} \; \begin{equation} \mathbb{G}_{l}(u) = \mathbb{H}_{l}(u) \; exp\left(-\dfrac{u^2}{2}\right), \; l=0,1,2,... \\ \end{equation} $ $ \begin{equation} \textbf{Hermite polynomial} \; \mathbb{H}_{l+1}(u) = 2u\mathbb{H}_{l}(u) - 2l\mathbb{H}_{l-1}(u) ,\:\text{where} \; \mathbb{H}_{0}(u)=1, \; and \; \mathbb{H}_{1}(u)=2u \end{equation} $

def hermite(x, order):
    if order == 0:
        return np.ones_like(x)
    elif order == 1:
        return 2 * x
    
    return 2 * x * hermite(x, order-1) - 2 * (order-1) * hermite(x, order-2)


def hermite_gauss(x, order):
    h = hermite(x, order)
    return h * np.exp(-x ** 2 / 2)

fig, ax = plt.subplots(1, 4, figsize=(32,5))

x = np.linspace(-6, 6, 500)

for i in range(4):
    hg = hermite_gauss(x, order=i)
    ax[i].plot(x, hg, color='red')
    ax[i].fill_between(x, hg, color='#D4D4D4')
    ax[i].grid(True)
    ax[i].set_title(f'HG with order {i}', fontsize=20)
    
plt.show()





# for i, _z in enumerate(z_list):
#     I = intensity(rho, _z, wavelength, z0, I0=1)
#     ax[i].imshow(I/I.max(), 
#                  cmap=cm, 
#                  extent=[x.min(), x.max(), y.min(), y.max()], aspect='equal')
#     ax[i].set_title(f"z={_z:.3f} mm", fontsize=20)

wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)


print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')


x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)

xx, yy = np.meshgrid(x, y)

rho = np.sqrt(xx**2 + yy**2)

Beam waist: 0.014192 mm
Divergence angle: 0.028385 rad
Depth of focus: 2.000000 mm

Intensity

The intensity of the Hermite-Gaussian beam $I_{l,m} = \left| U_{l,m} \right|^2$ is,

$$ \begin{equation} I_{l,m}(x,y,z) = \left|A_{l,m}\right|^2 \left[ \dfrac{W_{0}}{W(z)}\right]^2 \mathbb{G}_{l}^2\left[\dfrac{\sqrt{2}x}{W(z)}\right] \mathbb{G}_{m}^2\left[\dfrac{\sqrt{2}y}{W(z)}\right] \end{equation} $$

def hg_intensity(x, y, z, l, m, wavelength, z0, I0=1):
    W0 = np.sqrt(wavelength * z0 / np.pi)
    _W = W(z, wavelength, z0)
    
    hg_x = hermite_gauss(np.sqrt(2) * x / _W, l)
    hg_y = hermite_gauss(np.sqrt(2) * y / _W, m)
    
    return I0 * (W0 / _W * hg_x * hg_y) ** 2
    

fig, ax = plt.subplots(4, 4, figsize=(32,32))


for l in range(4):
    for m in range(4):
        hg = hg_intensity(xx, yy, 0, l, m, wavelength, z0)
        ax[m,l].imshow(hg, cmap=cm)
        ax[m,l].set_title(f"HG with order {(l,m)}", fontsize=30)
        ax[m,l].set_xticks([])
        ax[m,l].set_yticks([])

plt.show()

Laguerre-Gaussian beams

The complex amplitude of Laguerre-Gaussian beam has the form,

$$ \begin{equation} U_{l,m}(\rho,\phi,z) = A_{l,m} \left[ \dfrac{W_{0}}{W(z)}\right] \left( \dfrac{\rho}{W(z)}\right)^{l} \mathbb{L}^{l}_{m}\left(\dfrac{2\rho^{2}}{W^{2}(z)}\right) exp\left(-\dfrac{\rho^{2}}{W^{2}(z)}\right) exp\left[-jkz - jk\dfrac{\rho^{2}}{2R(z)} \mp{jl\phi}+ j(l+2m+1)\zeta(z)\right] \end{equation} $$

$ \textbf{Generalized Laguerre polynomial} \; \begin{equation} \mathbb{L}^{l}_{m+1}(u) = \dfrac{(2m+1+l-u)\mathbb{L}^{l}_{m}(u) - (m+l)\mathbb{L}^{l}_{m-1}(u)}{m+1} , \\ \:\text{where} \; \mathbb{L}^{l}_{0}(u) = 1 \; and \; \mathbb{L}^{l}_{1}(u) = 1 + l - u \end{equation} $

def laguerre(x, l, m):
    if m == 0:
        return np.ones_like(x)
    elif m == 1:
        return 1 + l - x
    else:
        return ((2*m-1+l-x)*laguerre(x, l, m-1) - (m-1-l)*laguerre(x, l, m-2)) / m

x = np.linspace(-2, 10, 500)

fig = plt.figure(figsize=(8,6))

ax = plt.subplot(111)


for l in range(4):
    for m in range(4):
        ax.plot(x, laguerre(x, l, m)) 
        
ax.grid(True)
ax.set_ylim(-5, 10)

plt.show()

x = np.linspace(-5, 20, 500)

fig = plt.figure(figsize=(8,6))
ax = plt.subplot(111)

for m in range(6):
    ax.plot(x, laguerre(x, 0, m)) 
        
ax.grid(True)
ax.set_ylim(-10, 20)
plt.show()

Intensity

The intensity of the Laguerre-Gaussian beam $I_{l,m} = \left| U_{l,m} \right|^2$ is,

$$ \begin{equation} I_{l,m}(\rho,\phi,z) = \left|A_{l,m}\right|^{2} \left[ \dfrac{W_{0}}{W(z)}\right]^{2} \left( \dfrac{\rho}{W(z)}\right)^{2l} {\mathbb{L}^{l}_{m}}^{2}\left(\dfrac{2\rho^{2}}{W^{2}(z)}\right) exp\left(-\dfrac{2\rho^{2}}{W^{2}(z)}\right) \end{equation} $$

def lg_intensity(rho, z, l, m, wavelength, z0, I0=1):
    W0 = np.sqrt(wavelength * z0 / np.pi)
    _W = W(z, wavelength, z0)
    
    lg = laguerre(2 * (rho / _W) ** 2, l, m)
    
    return I0 * (W0 / _W * lg) ** 2 * np.exp(-2 * (rho / _W) ** 2) * (rho / _W) ** (2*l) 

wavelength = 632.8e-6 # (mm)
z0 = 1;
W0 = np.sqrt(wavelength * z0 / np.pi)


print(f'Beam waist: {W0:.6f} mm')
print(f'Divergence angle: {2 * np.sqrt(wavelength / (z0 * np.pi)):.6f} rad')
print(f'Depth of focus: {2 * z0:.6f} mm')


x = np.linspace(-3*W0, 3*W0, 1000)
y = np.linspace(-3*W0, 3*W0, 1000)
z = np.linspace(-5*z0, 5*z0, 1000)

xx, yy = np.meshgrid(x, y)

rho = np.sqrt(xx**2 + yy**2) 

Beam waist: 0.014192 mm
Divergence angle: 0.028385 rad
Depth of focus: 2.000000 mm

fig, ax = plt.subplots(4, 4, figsize=(32,32))


for l in range(4):
    for m in range(4):
        lg = lg_intensity(rho, 0, l, m, wavelength, z0)
        ax[m,l].imshow(lg, cmap='jet')
        ax[m,l].set_title(f"LG with order {(l,m)}", fontsize=30)
        ax[m,l].set_xticks([])
        ax[m,l].set_yticks([])

plt.show()

Bessel beams (nondiffracting beams)

The complex envelope of Bessel beam has the form,

$$ \begin{equation} A(x,y) = A_{m} J_{m}(k_{T}\rho) e^{jm\phi}, \; \text{where} \; k_{T}^2 + \beta^2 = k^2 \end{equation} $$

$ \textbf{Bessel function of the first kind and $m^{th}$ order} \; \begin{equation} J_{m}(x) = \sum_{n=0}^{\infty}{\dfrac{(-1)^{n}}{n! \Gamma(n+m+1)} \left( \dfrac{x}{2}\right)^{2n+m}} \end{equation} $

Thus, the complex amplitude is,

$$ \begin{equation} U(r) = A_{m} J_{m}(k_{T}\rho) e^{jm\phi}, \; \text{where} \; k_{T}^2 + \beta^2 = k^2 \end{equation} $$

If $m=0$, the wave has complex amplitude, $ \begin{equation} U(r) = A_{0} J_{0}(k_{T}\rho) e^{-j \beta z} \end{equation} $ ,and the intensity distribution $I(\rho, \phi, z) = \left|A_{0}\right|^{2} J_{0}^{2}(k_{T}\rho)$

import plotly.graph_objects as go

fig= go.Figure(data=go.Isosurface(
    x=[0,0,0,0,1,1,1,1],
    y=[1,0,1,0,1,0,1,0],
    z=[1,1,0,0,1,1,0,0],
    value=[1,2,3,4,5,6,7,8],
    isomin=2,
    isomax=6,
))

fig.show()